library(sfcr)
library(tidyverse)

sim_eqs <- sfcr_set(
  TXs ~ TXd,
  YD ~ W * Ns - TXs,
  Cd ~ alpha1 * YD + alpha2 * Hh[-1],
  Hh ~ YD - Cd + Hh[-1],
  Ns ~ Nd,
  Nd ~ Y / W,
  Cs ~ Cd,
  Gs ~ Gd,
  Y ~ Cs + Gs,
  TXd ~ theta * W * Ns,
  Hs ~ Gd - TXd + Hs[-1]
)

sim_ext <- sfcr_set(
  Gd ~ 20,
  W ~ 1,
  alpha1 ~ 0.6,
  alpha2 ~ 0.4,
  theta ~ 0.2
)

system.time(sim <- sfcr_baseline(
  equations = sim_eqs, 
  external = sim_ext, 
  periods = 100, 
  hidden = c("Hh" = "Hs"),
  method = "Broyden")
  )
#>    user  system elapsed 
#>   0.172   0.028   0.199

sim %>%
  filter(period %in% c(1, 2, 3, 50)) %>%
  select(period, Gs, Y, TXd, YD, Hs) %>%
  t() %>%
  round(digits = 0)
#>        [,1] [,2] [,3] [,4]
#> period    1    2    3   50
#> Gs        0   20   20   20
#> Y         0   38   48  100
#> TXd       0    8   10   20
#> YD        0   31   38   80
#> Hs        0   12   23   80

shock1 <- sfcr_shock(
  variables = list(
    Gd ~ 25
  ),
  start = 5,
  end = 50
)

sim2 <- sfcr_scenario(sim, shock1, 50)

sim2_long <- sim2 %>%
  pivot_longer(cols = -period)

sim2_long %>%
  filter(name == "Y") %>%
  ggplot(aes(x = period, y = value)) +
  geom_line()

sim2_long %>%
  filter(name %in% c("YD", "Cd", "Hh")) %>%
  ggplot(aes(x = period, y = value)) +
  geom_line(aes(linetype = name, color = name))


sfcr_set_index(sim_eqs) %>%
  filter(lhs == "Cd")
#> # A tibble: 1 × 3
#>      id lhs   rhs                          
#>   <int> <chr> <chr>                        
#> 1     3 Cd    alpha1 * YD + alpha2 * Hh[-1]

simex_eqs <- sfcr_set(
  sim_eqs,
  Cd ~ alpha1 * YDE + alpha2 * Hh[-1],
  Hd ~ Hd[-1] + YDE - Cd,
  YDE ~ YD[-1],
  exclude = 3
)


simex <- sfcr_baseline(simex_eqs, sim_ext, 50, hidden = c("Hh" = "Hs"))


shock2 <- sfcr_shock(
  variables = sfcr_set(alpha1 ~ 0.7),
  start = 5,
  end = 50
)

simex2 <- sfcr_scenario(simex, shock2, 50)

simex2_long <- simex2 %>%
  pivot_longer(cols = -period)

simex2_long %>%
  filter(name %in% c("Cd", "YD", "Hh")) %>%
  ggplot(aes(x = period, y = value)) +
  geom_line(aes(linetype = name))

tfm_sim <- sfcr_matrix(
  columns = c("Households", "Firms", "Government"),
  codes = c("h", "f", "g"),
  c("Consumption", h = "-Cd", f = "+Cs"),
  c("Govt. Exp.", f = "+Gs", g = "-Gd"),
  c("Factor Income", h = "W * Ns", f = "-W * Ns"),
  c("Taxes", h = "-TXs", g = "+TXd"),
  c("Ch. Money", h = "-d(Hh)", g = "d(Hs)")
)

sfcr_validate(tfm_sim, sim, which = "tfm")
#> Water tight! The transactions-flow matrix is consistent with the simulated model.

sfcr_sankey(tfm_sim, sim)

sfcr_dag_cycles_plot(sim_eqs, size = 10)

sfcr_dag_blocks_plot(sim_eqs)